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ABSTRACT 
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Classical novae are powered by thermonuclear runaways that occur on the white dwarf component of close binary systems. During 
these violent stellar events, whose energy release is only exceeded by gamma-ray bursts and supernova explosions, about 10""^ - 
10~^ M© of material is ejected into the interstellar medium. Because of the high peak temperatures attained during the explosion, 
T'peak ~ (1 - 4) X 10^ K, the ejecta are enriched in nuclear-processed material relative to solar abundances, containing significant 
amounts of ^^C, ^^N, and ^^O and traces of other isotopes. The origin of these metal enhancements observed in the ejecta is not well- 
known and has puzzled theoreticians for about 40 years. In this paper, we present new 2-D simulations of mixing at the core-envelope 
interface. We show that Kelvin-Helmholtz instabilities can naturally lead to self-enrichment of the solar-like accreted envelopes with 
material from the outermost layers of the underlying white dwarf core, at levels that agree with observations. 

Key words. (Stars:) novae, cataclysmic variables - nuclear reactions, nucleosynthesis, abundances - convection - hydrodynamics - 
instabilities - turbulence 
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1. Introduction 

The assumption of spherical symmetry in classical nova mod- 
els (and in general, in stellar explosions) excludes an entire se- 
quence of events associated with the way that a thermonuclear 
runaway (hereafter, TNR) initiates (presumably as a point-like 
ignition) and propagates. The first study of localized TNRs on 
white dwarfs was carried out by Shara (1982) on the basis of 
semianalytical models. He suggested that heat transport was 
too inefficient to spread a localized TNR to the entire white 
dwarf surface, concluding that localized, volcanic-like TNRs 
were likely to occur. But his analysis, based only on radiative 
and conductive transport, ignored the major role played by con- 
vection on the lateral thermalization of a TNR. 

The importance of multidimensional eflTects for TNRs in 
thin stellar shells was revisited by Fryxell & Woosley (1982). 
In the framework of nova outbursts, the authors concluded that 
the most likely scenario involves TNRs propagated by small- 
scale turbulence. On the basis of dimensional analysis and 
flame theory, the authors derived the velocity of the deflagra- 
tion front spreading through the stellar surface, in the form 
Vdef ~ (^/?Vconv/Tbum)^^^, whcrc hp is the pressure scale height, 
Vconv the characteristic convective velocity, and Tbum the charac- 
teristic timescale for fuel burning. Typical values for nova out- 
bursts yield Vdef ~ 10^ cm s"^ (that is, the flame propagates 
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halfway throughout the stellar surface in about ~ 1.3 days). 
Shear-driven mixing induced by accretion of matter possessing 
angular momentum was also investigated by Kutter & Sparks 
(1987), but their numerical simulations failed to obtain a strong 
enough TNR to power a nova outburst (see Sparks & Kutter 
1987). 

The first multidimensional hydrodynamic calculations of 
this process were performed by Shankar, Arnett & Fryxell 
(1992) and Shankar & Arnett (1994). They evolved an accreting, 
1.25 Mq white dwarf with a 1-D hydro code that was mapped 
into a 2-D domain (a spherical-polar grid of 25x60 km). The 
explosive event was then followed with a 2-D version of the 
Eulerian code PROMETHEUS. A 12-isotope network, ranging 
from H to ^^F, was included to treat the energetics of the ex- 
plosion. Unfortunately, the subsonic nature of the problem, cou- 
pled with the use of an explicit code (with a timestep limited 
by the Courant-Friedrichs-Levy condition), posed severe limita- 
tions on the study, which had to be restricted to very extreme 
(rare) cases, characterized by huge temperature perturbations of 
about ~ 100-600%, in small regions at the base of the envelope. 
The total computed time was only about 1 second. The calcula- 
tions revealed that instantaneous, local temperature fluctuations 
cause Rayleigh-Taylor instabilities. Their rapid rise and subse- 
quent expansion (in a dynamical timescale) cools the hot mate- 
rial and halts the lateral spread of the burning front, suggesting 
that such local temperature fluctuations are not important in the 
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initiation or early stages of the TNR. The study, therefore, fa- 
vored the local volcanic-like TNRs proposed by Shara (1982). 

Glasner & Livne (1995) and Glasner, Livne & Truran (1997; 
GLT97) revisited these early attempts using 2-D simulations 
performed with the code VULCAN, din arbitrarily Lagrangian 
Eulerian (ALE) hydrocode capable of handling both explicit and 
implicit steps. As in Shankar et al. (1992), a slice of the star 
(0.1 TT^^^), in spherical-polar coordinates with reflecting bound- 
ary conditions, was adopted. The resolution near the envelope 
base was around 5x5 km. As before, the evolution of an ac- 
creting, 1 Mq CO white dwarf was initially followed using a 
1-D hydro code (to overcome the early, computationally chal- 
lenging phases of the TNR), and then mapped into a 2-D do- 
main as soon as the temperature at the envelope base reached 
Tb ~ 10^ K. As in the previous works, the 2-D runs relied on a 
12-isotope network. The simulations showed a good agreement 
with the gross picture described by 1-D models (for instance, 
the critical role played by the -unstable nuclei ^^N, ^"^'^^O, and 
^^F, in the ejection stage, and consequently, the presence of large 
amounts of ^^C, ^^N, and ^^O in the ejecta). However, some re- 
markable difl'erences were also identified. The TNR was initi- 
ated by an ensemble of irregular, localized eruptions at the en- 
velope base caused by buoyancy-driven temperature fluctuations 
indicating that combustion proceeds in a host of many localized 
flames - not as a thin front - each surviving only a few sec- 
onds. Nevertheless, these authors concluded that turbulent difl'u- 
sion efficiently dissipates any local burning around the core, so 
the fast stages of the TNR cannot be localized and the runaway 
must spread through the entire envelope. In contrast to 1-D mod- 
els, the core-envelope interface was convectively unstable, pro- 
viding a source for the metallicity enhancement of the envelope 
by means of a Kelvin-Helmholtz instability — resembling the 
convective overshooting proposed by Woosley (1986). Efficient 
dredge-up of CO material from the outermost white dwarf layers 
accounts for ~ 30% metal enrichment of the envelope (the ac- 
creted envelope was assumed to be solar-like, without any pre- 
enrichment), in agreement with the inferred metallicites in the 
ejecta from CO novae (Gehrz et al. 1998). Finally, larger con- 
vective eddies were observed, extending up to 2/3 of the enve- 
lope height with typical velocities Vconv ~ 10^ cm s"^ Despite 
these diff'erences, however, the expansion and progress of the 
TNR towards the outer envelope quickly became almost spheri- 
cally symmetric, although the initial burning process was not. 

The results of another set of 2-D simulations were pub- 
lished shortly afterward by Kercek, Hillebrandt & Truran 
(1998; KHT98), which aimed to confirm the general behav- 
iors reported by GLT97, in this case with a version of the 
Eulerian PROMETHEUS code. A similar domain (a box of 
about 1800x1 100 km) was adopted, but using a cartesian, plane- 
parallel geometry to allow the use of periodic boundary condi- 
tions. Two resolution simulations were performed, one with a 
coarser 5x5 km grid as in GLT97, and a second with a finer 
1x1 km grid. The calculations used the same initial model as 
GLT97, and produced qualitatively similar but somewhat less 
violent outbursts. In particular, they obtained longer TNRs with 
lower Tpeak and Vejec, caused by large diff'erences in the convec- 
tive flow patterns. Whereas GLT97 found that a few, large con- 
vective eddies dominated the flow, most of the early TNR was 
now governed by small, very stable eddies (with /^ax ~ 200 km), 
which led to more limited dredge-up and mixing episodes. The 
authors attributed these discrepancies to the difl'erent geometry 
and, more significantly, to the boundary conditions adopted in 
both simulations. 



The only 3-D nova simulation to date was performed by 
Kercek, Hillebrandt & Truran (1999), adopting a computational 
domain of 1 800x 1 800x 1000 km with a resolution of 8x8x8 km. 
It produced flow patterns that were dramatically difl'erent from 
those found in the 2-D simulations (much more erratic in the 3-D 
case), including mixing by turbulent motions occurring on very 
small scales (not fully resolved with the adopted resolution) and 
peak temperatures being achieved that were slightly lower than 
in the 2-D case (a consequence of the slower and more limited 
dredge-up of core material). The envelope attained a maximum 
velocity that was a factor -100 smaller than the escape velocity 
and, presumably, no mass ejection (except for a possible wind 
mass-loss phase). In view of these results, the authors concluded 
that CO mixing must take place prior to the TNR, in contrast to 
the main results of GLT9'fl. 

In summary, two independent studies, GLT97 and KHT98, 
based upon the same 1-D initial model, reached nearly oppo- 
site conclusions about the strength of the runaway and its capa- 
bility to power a fast nova. The origin of these difl'erences was 
carefully analyzed by Glasner, Livne & Truran (2005), who con- 
cluded that the early stages of the explosion, prior to the on- 
set of the TNR - when the evolution is almost quasi-static - 
are extremely sensitive to the outer boundary conditions (see 
e.g., Glasner, Livne & Truran (2007), for a 2-D nova simula- 
tion initiated when the temperature at the envelope base is only 
5 X 10^ K). Several outer boundary conditions were examined. 
The study showed that Lagrangian simulations, where the enve- 
lope is allowed to expand and mass is conserved, are consistent 
with spherically symmetric solutions. In contrast, in Eulerian 
schemes with a "free outflow" outer boundary condition — the 
choice adopted in KHT98 — , the outburst can be artificially 
quenched. 

In light of these conundrums, a reanalysis of the role of 
late mixing at the core-envelope interface during nova outbursts 
seems mandatory. To this end, we performed an independent 2-D 
simulation, identical to GLT97 and KHT98, with another multi- 
dimensional hydrodynamic code to investigate whether mixing 
can occur in an Eulerian framework with an appropriate choice 
of the outer boundary conditions. 

2. Models and input physics 

The 2-D simulation reported in this paper used FLASH, a paral- 
lelized, hydrodynamical, Eulerian code based on the piecewise 
parabolic interpolation of physical quantities for solving the hy- 
drodynamical equations and with an adaptive mesh refinement 
procedure. FLASH also uses a monotonicity constraint (rather 
than artificial viscosity) to control oscillations near discontinu- 
ities, a feature shared with the MUSCL scheme of van Leer 
(1979). For consistency with GLT97 and KHT98, the same ini- 
tial model was used. The model was computed by GLT97 on the 
basis of a 1-D, implicit hydro code, assuming accretion of solar 
composition matter (Z = 0.02) onto the surface of a 1 Mq CO 
white dwarf at a rate of 5 x 10"^ Mq yr"^. The accumulation 
of matter in degenerate conditions drives a temperature increase 
in the envelope, resulting in a superadiabatic temperature gradi- 



^ Other multidimensional studies (Rosner et al. 2001, Alexakis et al. 
2004a,b) focused on the role of shear instabilities in the stratified fluids 
that form nova envelopes. They concluded that mixing can result from 
the resonant interaction between large-scale shear flows in the accreted 
envelope and gravity waves at the interface between the envelope and 
the underlying white dwarf. However, to account for significant mixing, 
a very high shear (with a specific velocity profile) had to be assumed. 
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ent and eventually convective transport. The initial model corre- 
sponds to the time when the temperature at the innermost enve- 
lope zone is ^ 10^ K. At this stage, the mass of the accreted en- 
velope reaches 2 x 10"^ Mq. This radial profile has been mapped 
onto a 2-D cartesian grid of 800x800 km and is initially relaxed 
to guarantee hydrostatic equilibrium. The initial computational 
grid comprises 112 radial layers (including the outermost part 
of the CO core) and 512 lateral layers. Calculations rely on the 
adaptive mesh refinement with a minimum resolution 1.6x1.6 
km (simulations with a finer resolution will be presented in a 
forthcoming publication). 

A reduced nuclear reaction network was used to compute the 
energetics of the explosion: it consists of 13 isotopes (^H, "^He, 

12,13(-^ 13,14,15^^ 14,15,16,17q^ 17p _ Q^r^gj KHT98 

- supplemented with ^^F to include the important ^^0(p, y^^F 
reaction), linked through a net of 18 nuclear processes (mainly, 
p-captures and -decays). Reaction rates are taken from Angulo 
et al. (1999) and some more recent updates (see Jose, Hernanz, 
& Iliadis 2006, Jose & Shore 2008, and references therein). 

Periodic boundary conditions were adopted at both lateral 
sides, while hydrostatic boundary conditions are fixed at both 
the bottom (reflecting) and the top (outflow) The set of bound- 
ary conditions are similar to those implemented in GLT97 and 
KHT98, but note that the outer computational grid adopted in 
GLT97 is Lagrangian instead of Eulerian (to follow the late ex- 
pansion stages of the TNR). Finally, energy transport is included 
using an efl'ective thermal difl'usion coefficient that includes ra- 
diative and conductive opacities (Timmes 2000). 

3. Results 

In GLT97, significant numerical noise was present at the on- 
set of their calculations that produced temperature fluctuations 
of about 10-20%. We introduced (just at the initial time-step) a 
Gaussian temperature perturbation at the core-envelope interface 
of 5%. For comparison, the value in KHT98 was 1%. The size 
of the initial perturbation was 2 km, much smaller than the con- 
temporary depth of the accreted envelope (~ 800 km). The ini- 
tial perturbation produces fluctuations that move along the core- 
envelope interface during the first seconds of the simulations. 
These fluctuations, in turn, spawn Kelvin-Helmholtz vortices, 
which clearly show up about 200 s later (Fig. 1), and appear 
to initiate a turbulent cascade. Filaments and buoyant plumes 
are fully resolved in these simulations. At this stage, the fluid is 
characterized by a large Reynolds number, with a characteristic 
eddy length of 50 km, fluid velocities of 10^ - 10^ cm s"^ and 
a dynamic viscosity of 10^ R These Kelvin-Helmholtz (KH) in- 
stabilities transport unburnt CO-rich material from the outmost 
layers of the white dwarf core and inject it into the envelope. The 
characteristic eddy turnover time is Iconv/Vconv ~ 10 s. 

As the KH vortices grow in size, more CO-rich material is 
transferred into the envelope. Convection becomes more turbu- 
lent. The initially small convective eddies merge into huge shells 
(Fig. 2), as seen also in GLT97. At this stage, the nuclear energy 
generation rate reaches 10^^ erg g"^ s"^ while the characteris- 
tic burning timescale decreases to ~5 s. The convective filaments 
continue growing in size and progressively occupy the whole en- 
velope length. Although not resolved in these simulations, and 
in contrast to the 3-D case, the conservation of vorticity in 2- 
D forces the largest eddies to grow in an inverse vorticity cas- 
cade, while energy flows to the viscous scale with a distribution 

^ Technical details of how boundary conditions are implemented in 
Godunov-type codes can be found in Zingale et al. (2002). 



that deviates from the Kolmogorov spectrum (see e.g., Lesieur, 
Yaglom, & David 2000; Shore 2007). At this time, the tempera- 
ture at the envelope base reaches ~ 2 x 10^ K, at fluid velocities 
of 10^ cm s"^ (of the order of the escape velocity, characteristic 
of the dynamic phases of the explosion), and a nuclear energy 
generation rate of 10^^ erg g"^ s"^ . The convective turnover time 
is now ~ 5 s. The mean CNO abundance in the envelope has in- 
creased to 0.30, a value that agrees well with both the previous 
simulations by GLT97 and the mean metallicities inferred from 
observations of the ejecta in non-neon (CO) novae (see Jose & 
Shore 2008). At this stage, since the outer envelope layers had 
started to escape the computational (Eulerian) domain, simula- 
tions were stopped. 

Our 2-D simulations, in agreement with the results reported 
in GLT97, show that the progress and extension of the TNR 
throughout the envelope occurs with almost spherical symmetry, 
even though the structure of the ignition is not. This explains the 
success of 1-D models in reproducing the gross observational 
properties (light curves, velocities of the ejecta, nucleosynthe- 
sis) of nova explosions (Starrfield et al. 1998, 2009; Kovetz & 
Prialnik 1997; Yaron et al. 2005; Jose & Hernanz 1998). 

4. Discussion and Conclusions 

We have analyzed the possible self-enrichment of the solar- 
composition accreted envelope with material from the under- 
lying white dwarf during nova outbursts in a multidimensional 
framework. We have found that a shear flow at the core-envelope 
interface (which unlike the spherically symmetric case, does not 
behave like a rigid wall) drives mixing through KH instabilities. 
Large convective eddies develop close to the core-envelope inter- 
face, of a size comparable to the height of the envelope (similar 
to the pressure scale height in 1-D simulations), mixing CO-rich 
material from the outermost layers of the underlying white dwarf 
into the accreted envelope. The metallicity enrichment achieved 
in the envelope, Z ~ 0.30, is in agreement with observations of 
CO nova ejecta. Our 2-D simulations also show that even for a 
point-like TNR ignition, the expansion and progress of the run- 
away is almost spherically symmetric for nova conditions. We 
note that the adopted resolution as well as the size, intensity, and 
location of the initial perturbation have a very limited impact on 
the results, principally affecting the timescale for the onset of the 
KH instability but not the final, mean metallicity. Details will be 
extensively discussed in a forthcoming publication. Our results 
agree with earlier 2-D hydrodynamic simulations (GLT97) and 
solve the controversy raised by another 2-D study (KHT98) that 
questioned the efficiency of this mixing mechanism, and hence 
the corresponding strength of the runaway and its capability to 
power a fast nova outburst. 
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Fig. 1. Snapshots of the development of KH instabilities at ^ = 215 s (upper left panel), 235 s (upper right), 279 s (lower left), and 
498 s (lower right), shown in terms of ^^C mass fraction (in logarithmic scale). The injection of core material driven by the KH 
instabilities translates into a mass-averaged abundance of CNO-nuclei in the envelope of 0.079, 0.082, 0.089, and 0.17, respectively. 
The mean CNO abundance at the end of the simulations reaches 0.30, by mass. 




Fig. 2. Left and central panels: same as Fig. 1, but for the velocity fields at t = 279 s (left) and 498 s (right), superimposed on a plot 
of the nuclear energy generation rate (in ergs g"^ cm"^). Right panel: Time evolution of the overall nuclear burning rate. 
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